


# This is the empirical example for the paper "Subsampling Inference on Quantile Regression Processes" by V. Chernozhukov and I. Fernandez-Val
#NB: these programs are specialized to my empirical example

# Updated on 02/14/2014 to remove some commands and packages that have been superseded in R version 3.0.1

#Modified on 02/09/2021 to fit data used in SCLB 

############
rm(list=ls())

setwd("F:/Juan Munoz Dropbox/Juan Munoz/15. Papers UIUC/1. Mango Tree/MT_Juan/Juan's Folder/ReplicationFiles")
library(foreign)

table.rq.fn<- 
  function (formula, taus = c(0.25, 0.5, 0.75), method = "fn", data=data,
            ...) 
  {
    m <- length(taus)
    tab <- NULL
    for (i in 1:m) {
      fit <- rq(formula, taus[i], method = "fn", data=data)
      tab <- rbind(tab, coefficients(fit))
      
    }
    invisible(tab)
    
  }


table.rq.res<- 
  function (formula, taus = c(0.25, 0.5, 0.75), method = "fn", data=data,
            ...) 
  {
    m <- length(taus)
    tab <- NULL
    setab<- NULL
    for (i in 1:m) {
      fit <- rq(formula, taus[i], method = "fn", data=data)
      tab <- rbind(tab, coefficients(fit))
      setab <- rbind(setab, coef(summary(fit))[,2] )		
    }
    return(list(tab = tab, setab = setab))
    
  }



############## computes the inference processes for i-th subset of data ##############################

estimates<- function(i,ss, data, taus, formula){  
  
  datas<- data[ss[,i], ]
  if (det(t(as.matrix(datas)) %*% as.matrix(datas))>0.1)
  {;
    #sol<- t(table.rq(formula, taus=taus, method="br", data=datas)$a[,,"coefs"])	;
    sol<- table.rq.fn(formula=formula, taus=taus, method="fn", data=datas);
    solqr<-  sol[,2] ;
    solols<- lm(formula=formula,data=datas)$coef[-1]	;
    inference.process.shift<- solqr-solols[1]	; 
    inference.process.dom<-   (-(solqr > 0) + (solqr < 0)) * solqr ; 
    inference.process.noeffect<- solqr; 
    infer<- as.numeric(c(  	inference.process.shift,       inference.process.dom, inference.process.noeffect));
    write.table(t(as.matrix(infer)), file="SCLB.em.res_MT.txt", sep="  ", append=T,   col.names = F, na = "NA");
  };
  
}

############## subsampling function that computes the processes i \leq B bootstrap subsamples ############################

subsample<- function( data, nn=nn, b=b, B, formula) {
  ss<- matrix(sample(nn, b * B, replace = T), b, B);
  outpt<- lapply (1:B,  FUN=estimates, ss, data, taus, formula=formula)
  return(outpt)
}

##########################################################################################################################

# Part II:  Implementation

library(quantreg);
#	library(limma);
library(readstata13)
SCLB_data<- as.data.frame(read.dta13("DataSCLB.dta"));
keeps<- c("EL_EGRA_PCA_Index","MT_Program","CCT_Program", "BL_P1_w0", "miss_BL", "male", "Age")
SCLB_data<-SCLB_data[keeps]

#SCLB_data$group  <- factor(SCLB_data$Group_detailed)
#a           <- as.data.frame(model.matrix(~SCLB_data$group-1))
#colnames(a) <- (substring( names(a), 11, 24))
#SCLB_data <- cbind(SCLB_data, a)
#colnames(SCLB_data)[which(colnames(SCLB_data)=="group")] <- "group_FE"

nn<- dim(SCLB_data)[1];

#Cohort_FE <- paste("ct", 1:4, sep="_")
#Year_FE <- paste("yr", 1:4, sep="_")
#Group_FE <- paste("group", 1:45, sep="_")
#Group_detail_FE <- paste("g", 1:320, sep="_")

#formula <- formula(paste("EL_EGRA_PCA_Index ~ ", "MT_Program" ,"+CCT_Program",
                       #  "+BL_P1","+", paste(Cohort_FE, collapse=" + "),"+", paste(Year_FE, collapse=" + "),"+",
                        # paste(Group_FE, collapse=" + ")));

#formula <- formula(paste("EL_EGRA_PCA_Index ~ ", "MT_Program" ,"+CCT_Program",
#                         "+BL_P1","+", paste(Year_FE, collapse=" + "),"+",paste(Group_FE, collapse=" + ")));

#formula <- formula(paste("EL_EGRA_PCA_Index ~ ", "MT_Program" ,"+CCT_Program", "+BL_P1_w0", "+miss_BL", "+male", "+Age" ))

formula <- EL_EGRA_PCA_Index~MT_Program +CCT_Program +BL_P1_w0 +miss_BL


#ols results

z <- summary(lm(formula, data = SCLB_data));   ols <- z$coef;   vars <- names(z$coef[, 1]);  p <- length(z$coefficients[, 1]); Jn <- solve(z$cov.unscaled);


taus<- c(3:17)/20; Jt <- length(taus);
res <- table.rq.fn(formula, taus = taus, method = "fn", data = SCLB_data, hs=F ) ;

solqr<- res[,2]  ;
solols<- ols["MT_Program",1] ;

res.to.plot<- table.rq.res ( formula, taus=taus, method="fn", data = SCLB_data, hs=F    ) 


postscript("SCLB_QTE_MT.ps",horizontal=F, pointsize=15,width=7.0,height=7.0)

par(mfrow=c(1,1))

b<-  res.to.plot$tab[,"MT_Program"]   ;
b.p<- b + 1.69*res.to.plot$setab[,"MT_Program"] ; 
b.m<- b - 1.69*res.to.plot$setab[,"MT_Program"] ; 
plot( c(0,1 ),range(c( 0.12, b.m,b.p, -.12 )), xlim =c(0, 1), type="n",xlab="Quantile Index", ylab="Coefficient", main="Quantile Treatment Effect for Full Cost Program");
polygon(c(taus,rev(taus)),c(b.p,rev(b.m)),density=-3, border=T, col=3);
points(taus,b);
lines(taus,b);
abline(h=0);


dev.off()


inference.processN.shift <- solqr -solols;
inference.processN.noeffect <- solqr;
inference.processN.dom <- (-(solqr > 0) + (solqr < 0)) * solqr;


# subsampling stage to determine critical value

set.seed(runif(1)*100);
set.seed(88);

B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)



# CENTERED SUBSAMPLING

out<- read.table("SCLB.em.res_MT.txt");
B <- min(10000, dim(out)[1])
out<- out[1:B,-1]

# out<- matrix(unlist(output), B, 3*Jt, byrow=T)

out.shift<- out[,1:Jt] 
out.dom<-  out[,(Jt+1): (2*Jt) ] ;
out.noeffect<- out[,(2*Jt+1): (3*Jt) ] ;

varr.shift<- b*apply(out.shift, 2, var); 
varr.dom<-  b*apply(out.dom, 2, var); 
varr.noeffect<- b*apply(out.noeffect, 2, var);  

vb.shift<- sqrt( b*(out.shift -matrix(inference.processN.shift, B, Jt, byrow=T))^2/matrix(varr.shift, B, Jt, byrow=T) );
vb.noeffect<- sqrt( b*(out.noeffect -matrix(inference.processN.noeffect, B, Jt, byrow=T))^2/matrix(varr.noeffect, B, Jt, byrow=T));
vb.dom<- sqrt(b)* (out.dom -matrix(inference.processN.dom, B, Jt, byrow=T))/sqrt(matrix(varr.noeffect, B, Jt, byrow=T)) ;

vb.shiftN<- max(sqrt( nn*(inference.processN.shift)^2/varr.shift ));
vb.noeffectN<- max(sqrt( nn*(inference.processN.noeffect)^2 / varr.noeffect ));
vb.domN<- max(sqrt(nn)* (inference.processN.dom) /sqrt(varr.dom)) ;

max0<- function(y) { max(y, 0) }
vb.dom2<- apply(vb.dom, c(1,2),  max0)^2;

#  KSM.shift<- apply(vb.shift, 1, mean);  KS.shift<-  sqrt(nn*mean(inference.processN.shift^2/varr) );    
#  KSM.noeffect<- apply(vb.noeffect, 1, mean);  KS.noeffect<-  sqrt(nn*mean(inference.processN.noeffect^2/varr) );    
#  KSM.dom<- apply(vb.dom2, 1, mean);  KS.dom<-  sqrt(nn*mean(vN.dom2/varr) );    

KS.shift<- apply(vb.shift, 1, max);      
KS.noeffect<- apply(vb.noeffect, 1, max);      
KS.dom<- apply(vb.dom, 1, max);      

KSM.shift<- apply(vb.shift, 1, mean);      
KSM.noeffect<- apply(vb.noeffect, 1, mean);      
KSM.dom<- apply(vb.dom, 1, mean);      

# Unadjusted critical values

critical.shift<- quantile(KS.shift, .95); 
critical.noeffect<- quantile(KS.noeffect, .95); 
critical.dom<- quantile(KS.dom, .95); 



# record the KS, critical value, and the decision
reject.shift<- ifelse(vb.shiftN>critical.shift, "Reject", "Accept" )
reject.noeffect<- ifelse(vb.noeffectN>critical.noeffect, "Reject", "Accept" )
reject.dom<- ifelse(vb.domN>critical.dom, "Reject", "Accept" )


result.table <- rbind( 
  c(vb.shiftN,  critical.shift,  reject.shift  )
)   

result.table



write.dta(data.frame(result.table), file = "Output/AppendixTable7_FullProgram.dta")
